/*
 * EBFV1_assemblyRHS.cpp
 *
 *  Created on: 29/09/2012
 *      Author: rogsoares
 */


#include "EBFV1_elliptic.h"

namespace PRS{
	double EBFV1_elliptic::assemblyRHS(TMesh* pTMesh_data){
		setRHS_Dirichlet(pMatVec->G,pMatVec->RHS);
		setRHS_Neumann(pTMesh_data,pMatVec->RHS);
//		printVectorToFile(pMatVec->RHS,"RHS_G.txt"); STOP();
		int ndom = pTMesh_data->getNumDomains();
		for (int dom=0; dom<ndom; dom++){
			ierr = MatMatMult(pMatVec->E[dom],pMatVec->F[dom],MAT_INITIAL_MATRIX,1.0,&pMatVec->EF); CHKERRQ(ierr);
//			printMatrixToFile(pMatVec->EF,"Matrix_EF.txt");
			setRHS_Dirichlet(pMatVec->EF,pMatVec->RHS);
//			printVectorToFile(pMatVec->RHS,"RHS_EF.txt");
//			STOP();
		}
		int nrows_G = pMatVec->nrows;
		const int* rows_G = pMatVec->rows;
		int ncols_G = pMData->getNum_GF_Nodes();
		const int* cols_G = pMData->get_idxFreecols_ptr();
		ierr = MatGetSubMatrixRaw(pMatVec->EF,nrows_G,rows_G,ncols_G,cols_G,PETSC_DECIDE,MAT_INITIAL_MATRIX,&pMatVec->EF_freeNodes);CHKERRQ(ierr);
		return 0;
	}

	int EBFV1_elliptic::setRHS_Dirichlet(Mat A, Vec RHS){
		int numGP = pMData->getNum_GP_Nodes();	// get number of global prescribed (dirichlet) nodes
		int m,n;

		// get from matrix A a submatrix SUB_A. SUB_A columns are those columns from A which will fill the RHS vector for all free A matrix rows.
		Mat SUB_A;
		ierr = MatGetSubMatrixRaw(A,pMatVec->nrows,pMatVec->rows,numGP,pMData->get_idxn_ptr(),PETSC_DECIDE,MAT_INITIAL_MATRIX,&SUB_A);CHKERRQ(ierr);
//		printMatrixToFile(A,"Matrix_A.txt");
//		printMatrixToFile(SUB_A,"Matrix_SUB_A.txt");
//		cout << "pMatVec->rows:\n";
//		for(int p=0;p<numGP;p++){
//			cout << pMData->get_idxn_ptr()[p] << endl;
//		}
//		cout << "numGF = " << numGF << "\tnumGP = " << numGP << endl;
		//printVectorToFile(RHS,"RHS_G.txt");STOP();
		ierr = MatGetOwnershipRange(SUB_A,&m,&n);
		// loop over all local rows o submatrix SUB_A
		for (int row=m; row<n; row++){
			double sum = .0;
			MIter prescribedIter = pMData->dirichletBegin();
			for (int col=0; col<numGP; col++, prescribedIter++){
				double val = .0;
				ierr = MatGetValues(SUB_A,1,&row,1,&col,&val); CHKERRQ(ierr);
				sum += -val*prescribedIter->second;
			}
			//printf("sum[%d] = %f\n",row,sum);
			ierr = VecSetValue(RHS,row,sum,ADD_VALUES); CHKERRQ(ierr);
		}
//		VecView(RHS,PETSC_VIEWER_STDOUT_WORLD);
		ierr = VecAssemblyBegin(RHS); CHKERRQ(ierr);
		ierr = VecAssemblyEnd(RHS); CHKERRQ(ierr);
//		printVectorToFile(RHS,"RHS_G.txt"); STOP();
		ierr = MatDestroy(SUB_A); CHKERRQ(ierr);
		return 0;
	}

	int EBFV1_elliptic::setRHS_Neumann(TMesh* pTMesh_data, Vec RHS){
		int m,n,row,flag,ID;
		int numIW = pTMesh_data->getNumInjectionWells();
		for (int i=0; i<numIW; i++){
			flag = pTMesh_data->getInjectionWellFlag(i);
			// Do not include well flux on nodes with prescribed pressure
			if (pSimPar->isNodeFree(flag)){
				pTMesh_data->injectionWells_list[i].getSize(m,n);
				for (int j=0; j<m; j++){
					VertexWell *vwell = pTMesh_data->injectionWells_list[i].getValue(j,0);
					ID = vwell->vertex->ID_PETScOrdering - 1;
					row = pMData->FPArray(ID);
//					printf("%d\tQi[%d] = %.8E\n",flag,row,vwell->Qi);
					ierr = VecSetValue(RHS,row,vwell->Qi,ADD_VALUES); CHKERRQ(ierr);
				}
			}
		}
		int numPW = pTMesh_data->getNumProductionWells();
		for (int i=0; i<numPW; i++){
			flag = pTMesh_data->getProductionWellFlag(i);
			// Do not include well flux on nodes with prescribed pressure
			if (pSimPar->isNodeFree(flag)){
				pTMesh_data->productionWells_list[i].getSize(m,n);
				for (int j=0; j<m; j++){
					VertexWell *vwell = pTMesh_data->productionWells_list[i].getValue(j,0);
					flag = vwell->vertex->flag;
					ID = vwell->vertex->ID_PETScOrdering - 1;
					row = pMData->FPArray(ID);
//					printf("%d\tQi[%d] = %.8E\n",flag,row,vwell->Qi);
					ierr = VecSetValue(RHS,row,vwell->Qi,ADD_VALUES); CHKERRQ(ierr);
				}
			}
		}
		return 0;
	}
}

